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Many problems in robust control and motion planning can be reduced to either find a sound 
approximation of the solution space determined by a set of nonlinear inequalities, or to the "guar- 
anteed tuning problem" as defined by Jaulin and Walter, which amounts to finding a value for 
some tuning parameter such that a set of inequalities be verified for all the possible values of 
some perturbation vector. A classical approach to solve these problems, which satisfies the strong 
soundness requirement, involves some quantifier elimination procedure such as Collins' Cylindrical 
Algebraic Decomposition symbolic method. Sound numerical methods using interval arithmetic 
and local consistency enforcement to prune the search space are presented in this paper as much 
faster alternatives for both soundly solving systems of nonlinear inequalities, and addressing the 
guaranteed tuning problem whenever the perturbation vector has dimension one. The use of these 
methods in camera control is investigated, and experiments with the prototype of a declarative 
modeller to express camera motion using a cinematic language are reported and commented. 

Categories and Subject Descriptors: D.3.3 [Programming Languages]: Language Constructs 
and Features — Constraints; G.1.0 [Numerical Analysis]: General — Interval arithmetic; H.5.1 
[Information Interfaces and Presentation]: Multimedia Information Systems — Animations 

General Terms: Algorithms, Theory, Reliability 
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1. INTRODUCTION 

Designing electronic circuits [Ebers and JVIoU 1954], identifying the structure of com- 
plex molecules [Emiris and Mourrain 1999], or computing the quantity of chemical 
elements produced by some reaction [Meintjes and Morgan 1990] are all problems — 
among many others — that can be modelled by sets of real nonlinear equations and 
inequations. Reliably solving these systems (that is, delivering approximate solu- 
tions as close as possible to the true ones) with the limited set of real numbers 
representable on computers has generated a large amount of literature since the 
very dawn of computer science [Turing 1948; Rademacher 1948; Wilkinson 1963]. 
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Following Sunaga's and Moore's seminal works [Sunaga 1958; Moore 1966], inter- 
val analysis has been identified as a key tool for ensuring reliability and completeness 
(that is, the capability to retain all the solutions) in the solving process. 

Interval constraint solving [Older and Vellino 1990; Benhamou 1995] uses interval 
arithmetic in algorithms alternating propagation steps [Waltz 1975; Mackworth 
1977], which enforce some local consistency notion to tighten the variables' feasible 
domain, and search steps to isolate different solutions and overcome the weakness 
of the propagation steps. 

Interval constraint solvers such as clp(BNR) [Benhamou and Older 1997], ILOG 
Solver [Puget 1994], or Numerica [Van Hentenryck et al. 1997] have been shown to 
be efficient tools for solving some challenging nonlinear constraint systems [Puget 
and Van Hentenryck 1998; Granvillicrs and Benhamou 2001]). Relying on interval 
arithmetic, they guarantee completeness and isolate punctual solutions with an 
"arbitrary" accuracy. They take as input a constraint system and a Cartesian 
product of domains for the variables occurring in the constraints; their output is a 
set So of boxes approximating each solution contained in the input box. 

However, soundness (i.e., the property that all boxes returned contain nothing 
but true solution points) is not guaranteed while it is sometimes a strong require- 
ment. Consider, for instance, a civil engineering problem [Sam 1995] such as floor 
design where retaining non-solution points may lead to a physically infeasible struc- 
ture. As pointed out by Ward et al. [Ward et al. 1989] and Shary [Shary 1999], 
one may expect different properties from the boxes composing So depending on the 
problem at hand, namely: every element in any box is a solution, or there exists at 
least one solution in each box. Interval constraint solvers ensure only, at best, the 
second property. 

Interval constraint solving techniques also suffer from a severe limitation in that 
they can only handle constraint systems with discrete solutions. Yet, many ap- 
plications in robust control or error-bounded estimation [Jaulin and Walter 1993] 
are modelled by systems of nonlinear inequalities, whose solution set is usually not 
discrete. 

In addition, even more applications can be reduced to what Jaulin and Wal- 
ter [Jaulin and Walter 1996] call the "guaranteed tuning problem" , which amounts 
to finding the values for some tuning parameter such that a set of inequalities be 
verified for all the possible values of some perturbation vector^. More precisely: 

Given B, a box of feasible values for some tuning parameter vector 7 
and D, a box of feasible values for some perturbation vector tt, find the 
set S~f defined by: 

5^ = {7 g B I Vtt e D: f(7,7r) > 0} 

where f is a vector of nonlinear functions, and the inequality is to be 

taken com.ponentwise. 

Problems ranging from robust control [Abdallah et al. 1996] and camera con- 



-"^ Jaulin and Walter restrict the problem to finding only one value for the tuning parameter. The 
generalization adopted here is of interest to any application in which the user has to be given the 
choice of the solution to adopt, such as in the camera control application presented in this paper. 
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Fig. 1. Avoiding collisions with the arm of a robot 

trol [Druckcr and Zcltzcr 1994] to motion planning [Tarabanis 1990] can all be 
formulated as guaranteed tuning problems. Consider for example the following 
application: 

Example 1.1 A collision-free problem. A mobile robot arm is composed of three 

segments of respective lengths di, d2 and rfa, and three motor-controlled axes (see 
Figure 1). The trajectory of the robot's hand P{t) is therefore determined by three 
angle functions a\{t), 0:2 (i) and 0:3 (t). The problem is to find all points in the 
2D space that do not collide with the robot's hand, i.e. computing all the (.t, y) 
coordinates such that the distance between for all t, and (x, y) is greater than 
some given value d (here, the size of the robot's hand): 

Vi e [0, 1] : ^ {x - P^{t) f + {y - Py{t)Y ^ d 

where Px{t) and Py{t) represent the coordinates of P{t) at time defined by: 

r Px{t) = di sinai(t) + d2Sui{ai(t) + a2{t) - tt) + ^3 sin(ai(i) + a2{t) + az{t)) 
\ Py{t) = di cosai(t) + d2 cos[ai{t) + a2{t) - tt) + d3Cos(ai(t) + a2{t) + as{t)) 

Until recently, the soundness issue and the presence of quantifiers called for sym- 
bolic methods and quantifier elimination procedures such as Cylindrical Algebraic 
Decomposition [Collins 1975] (CAD). Unfortunately, these techniques are either too 
slow, or limited to polynomial constraints. 

The advent of interval analysis led to the devising of simple and sound algorithms 
to solve systems of inequalities [Jaulin and Walter 1993; Garloff and Graf 1999] 
and the guaranteed tuning problem (restricted to the finding of only one value, 
though) [Jaulin and Walter 1996]. By and large, many of these algorithms are but 
sophisticated interval extensions of simple search procedures, meaning that are they 
computationally expensive, even for small problems. 

In this paper, we present sound algorithms that draw upon efiicient complete 
interval constraint solving pruning methods to soundly solve systems of inequalities 
and the guaranteed tuning problem for the case of unidimensional perturbation 
vectors. This kind of problem is the one occurring in motion planning and camera 
control applications, where the only universally quantified variable is usually the 
time. 
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The outline of the paper is as foUows: in order to be reasonably self-content, the 
basics of interval analysis are presented in Section 2; related works on both solving 
systems of nonlinear inequalities and the guaranteed tuning problem arc presented 
in Section 3; their strong points and weaknesses regarding the applications targeted 
are also pointed out; interval constraint solving is introduced in Section 4 as a 
basis for the new soimd algorithms presented in Section 5; the modelling of a 
camera control problem in terms of these new sound interval constraint methods is 
then described in Section 6; the results with a prototype of a declarative modeller 
allowing a non-technician user to control the positioning of a camera by means of 
a cinematic language are commented in Section 7 and contrasted with the ones 
obtained by Jardillier and Languenou [Jardillier and Languenou 1998] on the same 
problems with a different approach; finally, Section 8 discusses directions for future 
researches. 

2. INTERVAL ANALYSIS 

In this section, we limit the presentation of interval analysis to the concepts needed 
in the sequel of the paper. The reader is referred to works by Moore, Hansen and 
others [Moore 1966; Alefeld and Herzberger 1983; Hansen 1992; Neumaier 1990] for 
a more complete presentation. 

The finite nature of computers implies that they can only represent and manipu- 
late a small subset of the real numbers represented in a floating-point format. Since 
the mid-eighties, most computers comply with the IEEE 754 standard [IEEE 1985] 
specifying the format of floating-point numbers. 

The set of floating-point numbers F is but a very small subset of real numbers. 
In addition, it is not closed for arithmetic operations, meaning that rounding of 
the results to representable floating-point numbers must usually take place. We 
say that an operation is correctly rounded when the value chosen, whenever the 
true result is not representable, is the closest floating-point number. Note that 
the IEEE 754 standard only requires the operators -h, — , x , ^ to be correctly 
rounded. The precision of the other operators is implementation dependent. See 
the paper by Lefevre et al. [Lefevre et al. 1998] for more information on this topic. 

Given a floating-point number a, let a+ (resp. a~) be the smallest float greater 
than a (resp. greatest float smaller than o). 

Rounding makes the reliable solving of systems of nonlinear equalities or in- 
equalities a challenging task, to which a large amount of papers and books has 
been devoted. 

One solution advocated by Moore [Moore 1966] to control the rounding errors 
is to use interval arithmetic, that is to replace reals by intervals containing them, 
whose bounds are representable numbers. For example: 

TT e [3.14,3.15] 

An interval I with representable bounds is called a floating-point interval. It is 
of the form: / = [a, 6] = {r € R | a ^ r ^ 6, with a,b €¥}. A non-empty interval 
I = [a, b] such that 6 < a"*" is said canonical. In the same way, a Cartesian product 
of intervals (or box) is said canonical whenever it is canonical in all its dimensions. 
In the sequel, vectors or Cartesian products are written in bold face. 
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Let I be the set of all floating-point intervals. Since we will only deal with this 
kind of intervals in the rest of the paper, we will refer to them as simply intervals. 

Operations and functions over reals are also replaced by an interval extension 
having the containment property. Given a real function, /: M" R and a box 
B e r , let /(B) = {/(r) | r e B}, Vf be the domain of /, and Vj = {B eT \ 
B C Vf}. 

Definition 2.1 Interval extension. Given a real function /: M" — » M, an interval 
extension f : I" — > I of / is an interval function verifying: 

Outer(/(B)) C F(B) VB e V) 

where Outer(p) = n{B G I" | p C B} for any real relation p C R". 

The extensions of some basic operators are defined as follows: 

[a, b] + [c, d]= [a + c,b + d] 

[a, h] — [c,d\= [a — d,h — c] 

[a, b] X [c, d] = [min(ac, ad, be, 6d), max(ac, ad, be, bd)] 

cxp([a. 6])= [cxp(a), exp(fe)] 

Note that in practice, the bounds computed have to be outward rounded in order 
to preserve the containment property. 

A particular interval extension, the natural interval extension is obtained by re- 
placing syntactically in a real function all the real constants by intervals containing 
them, all the real variables by interval variables, and all operators by their interval 
extension. 

Interval arithmetic is commutative, but it is neither associative (due to floating- 
point numbers wanting for this property themselves), nor distributive. It enjoys 
instead a sub- distributive property, namely: 

V/,J,iV'eI: I{J + K)CIJ + IK 

The sub-distributivity property has an important impact in that equivalent forms 
for a function over reals may not have equivalent natural extensions. 

In the rest of this paper, we will often consider the set of real numbers represented 

by a set of Cartesian product of intervals. Given Vi^S) the power set of any set <S, 
we then introduce the operator ^•]: 7^(1") K" defined as follows: 

VBi, . . . , VB„ e r : HBi, . . . , B„}5 = {r e R" | 3z € {1, . . . , n} s.t. r G B;} 

Lastly, given a box B = Ji x • • • x /„, an integer fc G {1, . . . , n}, and an interval 
J, let B/^<_j = 7i X • • • X 7fc_i X J X 7fe+i x • • • x 7„ be the box B where Ik has 
been replaced by J. 

3. RELATED WORK 

Computing a subpaving (Figure 2) of a real relation p C R" consists in partitionning 
R" into three sets {Ui,UoMu) of non-overlapping boxes with the properties: 

IU,1 C p 

lUo^np = (1) 

lUiUUu^Dp 
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Fig. 2. Subpaving of the relation p 



Given an n-ary constraint c, let pc Q K" be the underlying relation, that is: 



Given an n-ary constraint c and a box B G I", and assuming the existence of 
some procedure GlobSat with the following properties: 



it is easy to devise a systematic procedure to compute the subpaving of pc by 
alternating evaluation steps with GlobSat, and splitting steps whenever GlobSat 
returns "unknown." The precise algorithm is presented in Table I for a conjunction 
of constraints. 

The StoppingCriterion function appearing on Line 10 of Alg. Subpaving returns 
"true" or "false" depending whether the box given as an argument should still be 
considered for splitting. A typical instance of this method is testing for canonicity 
of the box. It is also possible to speed-up the computation by using another instance 
of StoppingCriterion that would avoid splitting boxes whose width is smaller than 
some £. 

The Splitj. function on Line 13 splits a box B into k non-overlapping subboxes 
whose union is equal to B. 

The 1+J operator on Line 14 applies on vectors of sets and performs their unions 
componentwise: 



The subpaving algorithm has been used by several authors to compute sound 
boxes for systems of inequalities. They differ by the way they implement the GlobSat 
method. 

In SIVIA, Jaulin and Walter [Jaulin and Walter 1993] use interval arithmetic 
containment properties: given a constraint c: /(xi, . . . , x„) ^ 0, and a box B, they 
evaluate the natural interval extension of / over B. If the right bound of the result 
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Table I. Subpaving algorithm for a conjunction of atomic con- 
straints ci A ■ ■ • A Cm 

1 Subpaving(in: {ci, . . .,Cm}, B e I"; out: {UiMoMu) 6 ^(1")^) 

2 begin 



3 sat ^ GlobSat(ci, B) A • • • A GlobSat(cm, B) 

4 switch sat in 

5 true; 

6 return ({B}, 0, 0) 

7 false: 

8 return (0, {B}, 0) 

9 unknown: 

10 if StoppingCriterion(B) then 

11 return (0, 0, {B}) 

12 else 

13 (Bi,...,Bk) ^Splitfc(B) 

k 

14 return |+) Subpaving({ci , . . . , Cm}, Bj) 

j=l 

15 endif 

16 end 



17 end 



is negative, GlobSat returns "true"; if the left bound is strictly positive GlobSat 
returns "false" ; otherwise it returns "unknown." SIVI A is able to process any kind 
of inequality constraints, be they linear or not. polynomial or not. The drawback 
of this approach is that for unstable / functions, the natural interval evaluation 
leads to large intervals that do not permit deciding whether the box B is included 
in pc or not. It is then necessary to split the box a lot. 

Jaulin and Walter [Jaulin and Walter 1996] have also devised an algorithm to 
solve the guaranteed tuning problem (restricted to the finding of only one value). 
It is based on SIVIA, and then, it suffers from the same drawbacks as SIVIA itself. 

The algorithm described by Kutsia and Schicho [Kutsia and Schicho 1999] im- 
plements another instance of Subpaving where GlobSat is obtained by testing some 
criterion using floating-point numbers of arbitrary precision. An important limita- 
tion is that their algorithm can only handle polynomial strict inequalities. 

Garloff and Graf [Garloff and Graf 1999] also restrict themselves to polynomial 
strict inequalities. They expand the polynomial inequalities into Bernstein polyno- 
mials: let I = (h, . . . be a multi-index (vector of non-negative integers), and 
= . . x^^ be a monomial. Given a polynomial p G R[xi, . . . , a;„], let <S be a 
set of multi-indices such that: p(x) = X]ie5 aiX"*". For any n-ary Cartesian product 
of domains D, one can express p in terms of Bernstein coefficients: 

p(x) = ^6i(D)i?^,i(x) 

where Sjv,i(x) is the Ith Bernstein polynomial of degree N. 

It is then well known that the following property does hold [Farouki and Rajan 
1987]: 

Va e D: min6i(D) < n(a) < max6i(D) (convex hull property) 
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As a consequence, it is possible to accurately bound the range of p(x) on every 
box, and then to know whether an inequality such as p{x) < does hold or not. 
By splitting the initial box and evaluating the corresponding Bernstein coefficients, 
one can isolate all the boxes verifying a set of inequality constraints. 

Sam-Haroud and Faltings [Haroud and Faltings 1994] do not rely on Alg. Sub- 
paving. Instead, they represent explicitly the subpaving of the relation p,. associated 
to each n-ary inequality constraint by a 2"-tree of boxes. The category of each box 
(in Hi, Uo, or is determined by finding the intersection of the curve with the 
edges of the box. Sets of constraints arc handled by performing intersections of 
the respective 2"-trees. In order to avoid a combinatorial explosion, n must be 
small. As a consequence, they have chosen to decompose the user's constraints into 
ternary constraints, so that they only manipulate octrees. The original algorithm 
is able to compute inner approximations (that is, subsets) of the solution space; on 
the other hand, it is not designed to handle the guaranteed tuning problem. 

The method presented by Collavizza et al. [Collavizza ct al. 1999] is strongly 
related to the one we present in the following, since they rely on usual interval 
constraint solving techniques to compute sound boxes for some constraint system. 
Starting from a seed that is known to belong to the solution space, they enlarge the 
domain of the variables around it in such a way that the new box computed is still 
included in the solution space. They do so by using local consistency techniques to 
find the points at which the truth value of the constraints change. Their algorithm 
is particularly well suited for the applications they target, viz. the enlargement of 
tolerances. It is however not designed to solve the guaranteed tuning problem. In 
addition, it is necessary to obtain a seed for each connected subset of the solution 
space, and to apply the algorithm on each seed if one is interested in computing 
several solutions (e.g. to ensure representativeness of the samples). 

In order to tighten a box B of variables' domains for a problem of the form 
Vu € Ik - ci A • • • A Cm, Jardillier and Languenou [Jardillier and Languenou 1998] 
compute an inner approximation by decomposing the initial domain Ik of v into 
canonical intervals ij..- ■ ■ ,1^, and testing whether ci A • • • A does hold for the 
boxes Ii X ■ ■ ■ X ll X ■ ■ ■ X In,- . . , /i x • • • x /|? x • • • x /„. These evaluations give 
results in a three- valued logic {true, false, unknown). Boxes labeled true contain 
only solutions, boxes labeled false contain no solution at all, and boxes labeled 
unknown are recursively split and re-tested until they may be asserted true or 
false, or canonicity is reached (see Figure 3). Retained boxes are those verifying: 

GlobSat(ci, Ji X - - - X ll X - ■ - X In) ^ true 
A 

A 

GlobSat(cTO, 7i X • • • X 7^ X • • • X 7„) = true 

The precise algorithm is described in Table II. In the initial call, the I parameter 
is equal to the left bound of v's domain. The Spllt^" procedure in Line 18 is identical 
to the Split^ procedure presented previously, except that it never splits the domain 
of the universally quantified variable v. 

Once again, Alg. JLA is but a sophisticated instance of Alg. Subpaving, which 

ACM Transactions on Computational Logic, Vol. V, No. N, February 2008. 



yj G{l,...,p}: { 



Interval Constraint Solving for Camera Control and Motion Planning • 9 



























V 



































































































a V b 



Fig. 3. JLA algorithm: solving Vj; G [a, b] : c{x,v) 

Table II. Evaluation-based propagation algorithm for Vi) G [v, tJ] : ci A ■ • • A Cm 

1 JLA(in: {ci, . . . ,c™}, B G I", a variable v,l&¥; out: {U^,Uo,Uu) G P{I"f) 
1 begin 



4 sat ^ GlobSat(ci, B) A ■ • ■ A GlobSat(cm, B) 

5 switch sat in 

6 true: 

7 if /+ = I) then 

8 return ({B}, 0, 0) 

9 else 

10 return JLA({ci , . . . , Cm}, B, v, Z+) 

11 endif 

12 false: 

13 return (0, {B}, 0) 

14 unknown: 

15 if StoppingCriterion (B) then 

16 return (0, 0, {B}) 

17 else 

18 (Bi,...,Bk) ^Spllt\"(B) 

k 

19 return |+| JLA({ci, . . . , Cm}, Bj, ti, i) 

j=i 

20 endif 

21 end 



22 end 



means that its efficiency strongly depends on the quahty of the GlobSat procedure. 
Like in SIVIA, JardiUier and Languenou use the natural interval extension of the 
constraints. As a consequence, their algorithm is computationally expensive in 
many cases, as soon as some moderate precision in the computation of the inner 
approximation is required. We refer the reader to Section 7 for precise figures. 

Due to lack of space, we will not present the Cylindrical Algebraic Decomposi- 
tion method for quantifier elimination here, and we refer the reader to the good 
introduction by Jirstrand [Jirstrand 1995]. 
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4. INTERVAL CONSTRAINT SOLVING 

Since finite representation of numbers by computers hinders the rehable solving of 
real equations and inequations, interval constraint solving relies on interval arith- 
metic [Moore 1966; Alefeld and Herzberger 1983] to compute verified approximate 
solutions of real constraint systems. 

The reader will find in various works [Older and VcUino 1990; Benhamou and 
Older 1997; Benhamou 1996; Van Hentenryck et al. 1997; Benhamou et al. 1994; 
Van Hentenryck et al. 1997; Benhamou et al. 1999] a thorough presentation of the 
interval constraint solving framework and of the associated state-of-the-art algo- 
rithms. In this section, we will describe only the elements that are essential to our 
purpose. In particular, we will restrict ourselves to the case of nonlinear inequality 
constraints only, that is constraints of the form f{xi, . . . ,Xn) ^ for some real 
function /. 

Proofs not given here may be found in the papers cited above. 
4.1 Local Consistencies 

Discarding all inconsistent values from a box of variables' domains is intractable 
when the constraints arc real ones (consider for instance the constraint c: sin(a;) = 
1, a; e [0, 2]). Consequently, weak consistencies have been devised, among which 
one may cite hull consistency [Benhamou 1995] and box consistency [Benhamou 
et al. 1994]. Both consistencies permit narrowing variables' domains to (hopefully) 
smaller domains, preserving all the solutions present in the input. Since they are 
used as a basis for the algorithms to be introduced in Section 5, they are both 
presented below. 

Discarding values of the variable domains for which c does not hold according 
to a given consistency notion is modelled by means of outer contracting operators, 
whose main properties are contractance, completeness, and monotonicity. 

4.1.1 Outer Contracting Operators. Depending on the considered consistency, 
one may define different contracting operators for a constraint. In this section, hull 
consistency and box consistency are first formally presented. Operators based on 
both consistencies are then given. 

Hull consistency is a strong "weak consistency" since a constraint c is said hull 
consistent w.r.t. a box whenever that box cannot be tightened without losing some 
solutions for c: 

Definition 4.1 Hull consistency. A real constraint c{xi, . . . ,Xn) is said hull- 
consistent w.r.t. a box B = /i x • • • x /„ if and only if: 

Vfc e {l,...,n}: 

Ik = Outer(/fe n {rfc e R I Vi G {1, . . . , n} \ {k} : 

3rj G Ij s.t. (ri,...,r„) G pc}) 

An operator enforcing hull consistency for a constraint c computes the smallest 
box containing the intersection of the input box and the relation Pf.: 

Definition 4.2 Outer-hull contracting operator. Let c be an n-ary constraint, pc 
its underlying relation, and B a box. An outer-hull contracting operator for c is a 
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function HCc: I" I" defined by: 

HCc(B) = Outer(Bn /9c) 

Proposition 4.3 Completeness of HC. Given a constraint c and a box B, 
the following relation does hold: (B n pc) C HCc(B). 

Operationally, hull consistency is enforced over constraints by decomposing them 
into conjunctions of primitive constraints. The drawback of such a method lies 
in the loss of domain tightening due to the introduction of new variables during 
the decomposition process. Box consistency has been introduced by Benhamou et 
al. [Benhamou et al. 1994] to overcome this problem: the operators enforcing box 
consistency consider constraints globally, without decomposing them. 

Definition 4.4 Box consistency. Let c be an n-ary real constraint, C an interval 
extension for c, and B = /i box. The constraint c is said box- consistent 

w.r.t. B if and only if: 

V/c G {!,..., n}: 

7fe = Outer(7fe n {r e M I C(/i, . . . , Outer({r}), /fe+i, . . . , /„)}) 

Intuitively, a constraint c is box-consistent w.r.t. a box B when each projection 
Ij,j G {l,...,n} of B is the smallest interval containing all the elements that 
cannot be distinguished from solutions of the univariate constraint obtained from 
C by replacing all the variables but Xj by their current domain due to the inherently 
limited precision of the computation with floating-point numbers. 

Box consistency is enforced over a constraint c{xi, . . . , Xn) as follows: n contract- 
ing operators (typically using an interval Newton method [Van Hentenryck et al. 
1997]) arc associated to the n univariate interval constraints obtained as described 
above. The fcth operator reduces the domain of Xk by computing the leftmost and 
rightmost canonical intervals J such that C(/i, . . . , Ik-i, J, Ik+i, ■ ■ ■ , In) does hold 
{leftmost and rightmost quasi-zeros). 

Using box consistency to narrow down the variable domains of a constraint leads 
to the notion of outer-box contracting operator: 

Definition 4.5 Outer-box contracting operator. Given an n-ary constraint c and 
a box B, an outer-box contracting operator BC3rc: I" — > I" for c is defined by: 

BC3rc(B) = max{B' | B' C B 

and c is box-consistent w.r.t. fc € {1, . . . , n} and B'} 

Proposition 4.6 Completeness of BC3r. Given a constraint c, the following 
relation does hold for any box B; 

(B n pe) C BC3rc(B) 

Consistencies and associated contracting operators considered so far are such 
that completeness is guaranteed (no solution lost during the narrowing process). 
Devising operators that ensure soundness of the results is the topic of the next 
section. 
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5. SOUND INTERVAL CONSTRAINT SOLVING 

We have seen in the previous section some operators that compute an outer ap- 
proximation of the relation p associated to some nonhnear constraint system. In 
this section, we will present operators that compute an inner approximation of p, 
that is, a subset of the solution space of the corresponding constraint system. As 
a consequence, the boxes computed by these operators will have the property that 
any point inside of them is a true solution to the problem, thus ensuring soundness 
of the results given to the user. 

Several definitions for the inner approximation of a real relation exist in the 
literature, depending on the intended application. Given an n-ary relation p, one 
may single out at least the two following definitions for an inner approximation 
Inner(p) of p: 

(1) Def. A. lnner(/3) = Bi, where Bi S {B S I" | B C p} (an inner-approximation 
is any box included in the relation) [Markov 1995; Armengol et al. 1998]; 

(2) Def. B. Inner(p) = Bi, where Bi G {B = /i x • • • x J„ | B C p A Vj e 
{1, . . . ,n},V/j 3 Ij : B/^<_j( C p I'- = Ij} (an inner- approximation is a 
box included in the relation that cannot be extended in any direction without 
containing non-solution points) [Shary 1999]. 

Definitions A and B imply that disconnected relations are only very partially 
represented by one box only, a drawback that is avoided with the following stronger 
definition, which will be used in this paper: 

Definition 5.1 Inner approximation operator. Given an n-ary real relation p C 
M", the inner approximation operator Inner: M" M" is defined by: 

lnner(/9) = {r e M" | Outer({r}) C p} 

The inner-approximation contains all the elements whose enclosing box is in- 
cluded in the relation. The Inner operator enjoys the following properties: 

Proposition 5.2 Properties of the Inner operator. The Inner operator is 

contracting, monotonic, idempotent, suh- distributive w.r.t. the union of subsets of 
K", and distributive w.r.t. the intersection of subsets o/K". 

Proof. Given A, B, and C three subsets o/M", with AC B, let r be an element 
o/R" and D = Outer({r}). 

Contractance. Given any real r e Inner(p), we have by definition of Inner that 
Outer({r}) C p. By definition of Outer, we have that r e Outer({r}), which leads 
to r e p. Gonsequently, Inner(p) C p; 

Monotonicity. For any r, r e lnner(>t) implies T) C A, by definition of Inner; 
since A Q B , we have D C B, and finally, by definition of Inner, r S Inner(B). 
Gonsequently, AC B Inner(^) C lnner(,B); 

Idempotence. We first prove lnner(lnner(^)) C Inner(^).- 

Given r G lnner(lnner(yl)), we have D C Inner(^), by definition of Inner. By 
contractance of Inner, Inner(^) C A. We then have D C and then r e Inner(^), 
by definition of Inner again. 

Now, we prove lnner(lnner(>l)) D lnner(>l); 
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Given r G Inner(^), we have T) (~ A hy definition of Inner. Therefore, any 
element in D is in Inner(^). As a consequence, D C Inner(^). Definition of Inner 
then leads to r d lnner(lnner(^)); 

Sub-distributivity w.r.t. union. We prove that: 

Inner(e) U Inner(C) C lnner(S U C) 

We have: 

Inner(S) U Inner(C) = {r e M" | Outer({r}) C B} U {r e M" | Outer({r}) C C} 
= {r e M" I Outer({r}) CBV Outer({r}) C C} 
C {r e M" I Outer({r}) CBUC} 
C Inner(^UC) 

Distributivity w.r.t. intersection. We prove that: 

lnner(6) n Inner(C) = Inner(i3 n C) 

We have: 

Inner(enC) = {r G M" | Outer({r}) C {BOC)} 

= {r e I Outer({r}) C B A Outer({r}) C C} 

= {r e M" I Outer({r}) C B} n {r e M" | Outer({r}) C C} 

= Inner(B) n Inner(C) 

□ 

5.1 Inner Contracting Operators 

The narrowing of variable domains occurring in a constraint is done in the same 
way as in the outer-approximation case: an inner contracting operator associated to 
each constraint discards from the initial box all the inconsistent values along with 
consistent values that cannot be enclosed in a canonical computer-representable 
box. The result is a set of boxes. 

Definition 5.3 Inner contracting operator. Let c be an n-ary constraint. An in- 
ner contracting operator for c is a function ICc : I" ^ 7^(1") verifying: 

VB e r : ^ICc(B)5 C lnner(B n pc) 

Proposition 5.4 Soundness of IC. Given a constraint c and an inner-con- 
tracting operator ICc for c, we have: 

VBel": ^ICc(B)5c (Bnpc) 

Proof. Immediate consequence of Inner and Outer definitions. □ 

Remark 5.5. Given a constraint c, an inner contracting operator ICc for c, an 
outer contracting operator OCc for c, and a box B, the following relations, deriving 
from Prop. 4.3 and Prop. 5.4, hold: 

^ICe(B)^ C (B n pe) C OCe(B) 
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Inner-contracting operators with stronger properties (computation of the greatest 
representable set included in a relation) may be defined, provided some strong 
assiimption (namely the ability to determine for any canonical box whether it is 
included or not in the relation), is fulfilled. These operators are optimal in the 
sense defined below. 

Definition 5.6 Optimal inner contracting operator. Let c be an n-ary constraint, 

and ICc an inner-contracting operator for c. The operator ICg is said optimal if and 
only if the following relation docs hold for any box B: 

^IC,(B)] = lnner(Bnp,) 

Let us consider the following naive algorithm as an implementation of an optimal 
inner contracting operator: given a constraint c, and a box B of domains for the 
variables occurring in c, let us split B in all its canonical subboxcs. We can apply 
the outer- hull contracting operator HCc enforcing hull consistency on each of these 
canonical subboxes. Consider the possible cases for one of the canonical subboxes 
D: 

o D n Pc = 0- By definition of hull consistency, we must have HCc(D) = 0; 
o D (1 pc ^ 0- There are again several cases (considered exclusively, from top to 
bottom): 

• T) C p^: By completeness of hull consistency-based operators, we must have 
HCc(D) = D. Now, if we consider the negation c of c, we must have HCc(D) = 
0, by definition of hull consistency, since D does not contain any element of 

Pc, 

■ D n /Oc ^ 0. The box D contains both some elements of Pc and Conse- 
quently, by completeness, HCc(D) and HCc(D) must return non-empty boxes 
included in D (a canonical box may still be tightened by shrinking any of its 
non-punctual dimensions to one of the bounds). 

We then have a procedure to determine for each canonical box D whether it is 
part of lnner(/9c) or not, namely: 

o HCc(D) = =^ D 2 Inner(pc); 
o HCc(D) = D C Inner(pc); 
<> (HCc(D) C D A HCc(D) CD) ^ D ^ Inner(pe). 

By contractance and completeness of HCc, there are no other possible outcomes 
(note that all cases are considered mutually exclusive for any non-empty canonical 
interval). 

Now, the preceding rules do not allow to implement an optimal inner contracting 
operator in practice for the following reasons (disregarding the fact that it would 
be grossly inefficient anyway): 

— We have restricted our presentation of consistency operators to closed interval 
arithmetic only. As a conscriucncc, for any constraint of the form /(xi, . . . , a;„) ^ 
0, we cannot consider its negation /(xi, . . . , a:„) > 0, which would require open 
interval arithmetic as well. Instead, we will define its negation as the constraint 
f{xi, . . . , Xn) > 0. As a result, the property "HCc(D) ^ D ^ Inner(pc)" 
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Table III. Inner contracting operator for an 
atomic constraint c based on the outer contracting 
operator F 

1 IC01^(in: B € I"; out: {UiMoMu) G P(n")3) 

2 begin 

3 B'^rc(B) 

4 Wo ^ B S B' 

5 B" ^ r^(B') 

6 Wi ^ B' H B" 

7 if StoppingCriterion(B") then 

8 return (Wi, Wo, {B"}) 

9 else 

(Bi,...,Bk)^Split^(B") 

11 return {Ui,Uo, 0) tti ^ (+) IC0l5^(Bj)^ 

12 endif ^ 

13 end 



docs not hold any more. Consider for example the constraint c: x ^ 0.5 with x € 
[0.5,0.5+]. We have HCc([0.5, 0.5+]) = [0.5,0.5+]; considering c as the constraint 
X < 0.5, we obtain: HCc([0.5, 0.5+]) = [0.5,0.5]. Applying the rules above, we 
would conclude that the interval [0.5,0.5+] does not belong to Inner(pc), though 
it does. It is however important to note that the "relaxed" definition for the 
negation we arc forced to adopt has no effect on the soundness of the algorithms 
to be described. On the other hand, it will preclude us from devising optimal 
inner contractors; 

— Even if wc were using open interval arithmetic, we still would have to face the 
fact that we are usually not able to implement operators enforcing hull consis- 
tency on the constraints given by the user but on primitives obtained after the 
decomposition process; 

— Lastly, wc have seen in Section 2 that the correct rounding of most floating-point 
arithmetic operators is usually not guaranteed, which precludes us from imple- 
menting contracting operators that enforce hull consistency, even for primitive 
constraints. 

Nevertheless, we can still use this principle to implement non-optimal inner con- 
tracting operators based on any kind of outer contracting operator. Table III 
presents such an implementation. The operator ICOl is an inner contracting oper- 
ator for a constraint c parameterized by an outer contracting operator T for c. 

The operator S occurring on lines 4 and 6 returns the set difference between two 
boxes as a set of boxes, that is: 

VBi,B2 G I": Bi SB2 C V{r) 

with ^Bi S Baj = {r e Bi | Outer({r}) n Ba = 0} 

The result is not uniquely defined since there are many ways to represent a set by 
unions of boxes. This is however irrelevant to our purpose. 

Proposition 5.7 Soundness of ICOl. Given an n-ary constraint c, a box's e 
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Table IV. Inner contracting operator for a constraint "iv & J : c 

1 IC02e (in: B G I"; out: {Ui,UoMu) € ^(I")^) 

2 begin 



3 B' ^ rc{B) 

4 if Domg/(^') = Dom-B{v) then 

5 iYo B H B' 
B" ^ r^(B') 

7 Wf ^ B'Domg,(i;)^J S B"Dorn3„(-u)^J 

8 if StoppingCriterion(B") then 

9 return (Wi, Wo, {B"}) 
10 else 

(Bi,...,Bk) ^-Split^^CB") ^ 

12 {UiMoMu) ^ {UiMo,0) td [9 "^°2c'^''"''(Bj) 

13 return {Ui,Uo,Uu) 

14 endif 

15 else 

16 return (0, {B}, 0) 

17 endif 



18 end 



I", an outer contracting operator T for c, and {]Ai,UoMu) = ICOlJj^(B), we have: 

( VDeUo-. D n pc = 
\VDeWi: DCp, 

Proof. Consider Line 3 of ICOl. By completeness of T , B n = B' n Pc- 
Consequently, ^B S B'jflpc = 0. Since the setUg returned eventually is the union 
of all the sets computed on Line 4, we have that Uof\Pc = ^- We can use the dual 
reasoning to prove that VD e ZY, : D C p^- □ 

5.2 Solving the Unidimensional Guaranteed Tuning Problem 

We now consider the Unidimensional Guaranteed Tuning Problem (UGTP): given 
a set of nonlinear inequalities ci,...,Cm on n variables {xi, . . . ,Xn-i,v}, a box 
Bel", and a domain J, compute the set I defined by: 

J= {re B I Vw e J: ci A--- Ac„} (2) 

If we restrict ourselves to the case m = 1, it is easy to devise an algorithm 
implementing an inner contracting operator computing a subset of I along the 
same lines as ICOl: given the constraint Vv G J: c, use ICOl to compute an inner 
approximation of pc while retaining only the boxes for which the domain of v is 
equal to J. The modified algorithm is presented in Table IV. The algorithm IC02 
is parameterized by an outer contracting operator F for c, by the variable v that is 
universally quantified, and by the domain J on which v must range. 

The operator DomB(i') occurring on Line 4 returns the interval in box B that 
represents the domain of v. Note that after applying Pc on Line 3, we check that 
the domain of v has not been tightened. Otherwise, there would exist some value 
in the domain of v such that c does not hold. The box B would then have to be 
discarded. The domain of v may however be tightened on Line 6 from some domain 
Ja to some domain Ji, by the outer contracting operator for c. This corresponds to 
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having proved that Vxi G \ /(' . . . Vx„_i G In-i \ I'n-i^v G Ja\Jb- c. It is then 
sufficient to prove Vu e Jfe : c in the next steps of the algorithm. 

In the following, given an n-ary constraint c, an interval J, and a variable v, let 
us note p^v the n-ary relation-^ associated to the constraint Vz; G J : c. 

We first state a simple lemma, which will be used in proving the proposition to 
follow. 

Lemma 5.8. Given n G N, and two boxes B = 7i x • • • x 7„ and D = Ji x • • • x J„ 
such that BCD, the following does hold: 

n 

D\B = |JJiX---x Jfe_i X (Jfe \ Ik) X Jfe+i X • • • X J„ 

k=l 

In addition: 

n 

S Bj = [J Ji X • • • X Jk-l X IJk S 45 X Jk+l X ■■■ X Jn 
k=l 

Proof. Immediate. □ 

Proposition 5.9 Soundness of IC02. Let c{xi, . . . ,Xn-i,v) be an n-ary con- 
straint, B = Ji X • • • X In-i X J a box, T an outer contracting operator for c, and 
{JAi,UoMu) = IC02j^'^"''^'' (B). The following properties do hold: 

{ yz) eUoi D n = 

Proof. Termination of Alg. IC02 depends on a reasonable choice for the func- 
tions StoppingCriterion() and Split^"(). More specifically, termination is ensured 
whenever Split^"() creates boxes Bi, . . . , B^ that are all strictly smaller than B", 
and when StoppingCriterion() checks for the size of all dimensions of but the 
one of V being smaller than some threshold that is attainable with the floating-point 
format at hand. 

In order to prove the proposition, we have to show that the properties given by 
Eq. (3) always hold on lines 9 and 16 in Table IV. Provided that boxes put into Uq 
andUi on lines 5 and 7, respectively, verify Eq. (3), it is not necessary to investigate 
what is returned on line 13 since it is the mere union of the sets previously computed. 

► We first consider the case when DomB'(") 7^ DomB(^')■ Let B = /i x • • • x /„. 
By hypothesis, DomB(v) = In = J during the first call of IC02. By completeness 
ofTc, i/DomB'(f) 5: DomB(w), there exists a value in DomB(w) \ DomB'(w) such 
that c does not hold when the free variables take their values in Ii x ■ ■ ■ x In-i- 
Consequently, B n = 0- The property does hold for all the subsequent calls of 
IC02 since the corresponding domains for v are always included in J (by contrac- 
tance of the narrowing operators used). Since C p^v , we have proved that the 
triplet returned on Line 16 verifies Eq. (3) for any call of IC02. 

► We now consider the case when DomB'(f) = DomB(i'). 

■^For convenience, the (n— l)-ary relation associated to Vd S J: c is extended to an n-ary relation 
by using the domain J for the extra dimension. 
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Let us prove that for any call I of IC02, the assignments performed on lines 5 
and 7 are such that: 

This suffices to prove the proposition since what is returned for the case DorriB' (v) = 
DomB(i') is either Uo^ andlAf \ or the union, componentwise, ofUo^ anduf^ with 
sets of boxes ui™'^ and U^"^^ (m > I) computed in subsequent calls, and since we 
have already proved that what is returned whenever DorriB' (v) ^ Dome (v) verifies 
the property. 

> We first consider Line 5: by completeness ofTc, S B'j C p^. Hence, lUo^ C 
Pc, and then JjUo^ H p^* = <ZS. 

> We now consider Line 7: we will first prove that for any call I of IC02; 

mjf^ .sB^'Hcpc (4) 

• This is trivially true for the first call (1 = 1) since, by hypothesis, DomB(i)(w) = J . 

• Assume that Eq. (4) does hold for the l-th call of IC02. Since B' C B and 
Domg,(i) (t^) = Domgd) (t^), we have: 

Kl,„,(.wSB'(')5Cp, (5) 
By completeness of pc, we have: 

^B"Dom3„(o (,;)^Dom^,(o (v) ^ ^"^'^^ C (6) 

By contractance ofT-c, B"''^ C B'^'\ From equations (6) and (5), we deduce: 

^B"Ll^„,,(„)^,SB"(')^Cp, (7) 

The box B"*-'-* is split on Line 11 in k boxes Bi^'^ (i & {1, . . . ,k}), with Domg.co (v) = 
DomB//(i)(u) (no splitting on the domain ofv). 
Consequently, we obtain from Eq. (7); 

Vie{l,...,fc}: ^BiW^^^^^^(„)^^sBiW5Cp, 

From Line 12, it follows that Eq. (4) does hold for I -\- 1, . . . ,1 + k whenever it does 
hold for I. 

Having proved that Eq. (4) does hold for any call I, we reconsider Line 7: by 
completeness ofVc, it comes: 

l^' s B"5 c Pc 

Hence: 

IB' S B" DomB„(j))^j5 C Pc (8) 

since. B"qqi^^„(„)^j D B". 

From Eq. (5) and (8), we deduce: 

^B'DomB,(u)^J S B"Don,B//(v)^J^ ^ Pc (9) 
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Table V. Inner propagation algorithm for a set of constraints S based 
on the inner contracting operator T 



1 IPA^(in: a set of constraints 5,B 6 I"; out: {Ui,UoMu) 6 ^(1")^) 



2 begin 

3 r^{B} 

4 (UoMu) {!3,<3) 



5 foreach c G S do 

6 Ui ^ 

7 foreach D g T do 



(UiMoMu) ^ {UiMoMu) w Tc(D) 



9 end 

10 T <—Ui 

11 end 



12 return {Ui,Uo,Uu) 
12 end 



We also know from Lemma 5.8 that all the boxes in ^'Dom^,{v)^j ^^"Dorr\^„{v)^j 
have J as their last dimension. From this and Eq. (9), we obtain: 



That is, lUi^ C p^v for any call of IC02. □ 

Handling the Unidimcnsional Guaranteed Tuning Problem described by Eq. (2) 
is done by the propagation algorithm I PA presented in Table V as follows: each 
constraint of the system is considered in turn together with the sets of elements 
verifying all the constraints already considered; the main point concerning I PA lies 
in that each constraint needs only be taken into account once, since after having been 
considered for the first time, the elements remaining in the variable domains are all 
sohitions of the constraint. As a consequence, narrowing some domain later does not 
require additional work. Alg. I PA is parameterized by an inner contracting operator 
T. Solving the UGTP can be done by instantiating IPA with IC02. Using ICOl 
instead loads to an algorithm computing the inner approximation of the relation 
associated to the conjunction of constraints in S. 

Proposition 5.10 Soundness of Alg. IPA. Given a set of constraints S = 
{ci, . . . ,Cm}, a box B e I", an inner contracting operator T, a variable v, an 
interval J, a set of outer contracting operators {Fi, . . . ,T„i} for, respectively, ci, 
Cm, and {Ui,UoMu) = IPA^(iS, B), let us note ps, the relation associated 
to the constraint ci A ■ ■ ■ A Cm and p^v the relation associated to the constraint 
Vv e J: Ci A • • • A Cm. We have: 



Proof. For each constraint c, we consider only the boxes that were put in Ui by 
thepreceding constraint (Lines 7 and 10). By soundness of ICOl (resp. \C02), when 
considering constraint Cj, Ui then contains on line 10 only the boxes D verifying 




D n PS = 

DCp5 

D n = 

D C p^v 



D C 



Pci A • • • A D C pc^ (resp. D C p^v A • • • A D C p^y/ 
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Once again, by soundness of ICOl (resp. \C02), when considering a constraint 
Cj on Line 5, Uo contains only boxes D such that there was a constraint Ck with 
k < j for which T) H pc^ =0 (resp. D n p^v = 0). Consequently, D n (pci H • • • fl 

Pc. n • • • A pe,) = 0) (resp. D n (p,v n • • • n n • • • a p,v) = 0). □ 

Note that it is straightforward to modify I PA in order to be able to solve con- 
straints of the form V?; e J^^^ : Ci A • • • A Vt; G J^™^ : Cm by replacing the set of 
constraints S by a set of pairs {ci,J^'^^) and by initializing accordingly the box 
passed to Tc on Line 8. Since we have not encountered applications requiring such 
extension, we will not consider it any more in the following. 

6. CAMERA CONTROL AND THE VIRTUAL CAMERAMAN PROBLEM 

Camera control is of interest to many fields, from computer graphics (visualization 
in virtual environments [Blinn 1988]) to robotics (sensor planning [Abrams and 
Allen 1997]), and cinematography [Davenport et al. 1991]. Whatever the activity, 
the objective is always to provide the user with an adequate view of some points of 
interest in a scene for a predefined duration. In sensor planning [Tarabanis 1990] 
for instance, one is interested in positioning a camera over a robot in order to be 
able to monitor its work whatever the position of its manipulating hand may be. 
The targeted applications are mostly real-time ones, which implies seeking for only 
one solution through some optimization process. The modelling adopted is usually 
very close to the camera representation (involving for instance, the direct control of 
the camera parameters through the use of sliders), thereby impeding inexperienced 
users from predicting the exact behaviour of the controlled device. 

Yet, some authors [Glcicher and Witkin 1992; Drucker et al. 1992] from computer 
graphics and cinematography fields have worked on determining camera parameters 
from given properties of a desired scene. Once more, most of these works are 
concerned with the computation of only one solution by means of an optimization 
criterion. 

Among these works, one may single out the attempt by Gleicher, which permits 
using some constraints including the position of a three-dimensional point on the 
screen, and the orientation of two points along with their distance on the projected 
image. Higher level of control is obtained through the ability to bound a point 
within a region of the image or to bound the size of an object. Time derivatives 
of the camera parameters are computed in order to satisfy user-defined controls. 
The method is devoted to maintaining user-defined constraints while manipulating 
camera parameters. Camera motion in an animated scene is not treated. 

Another approach [Jardillier and Languenou 1998] — inspired by Snyder's sem- 
inal work [Snyder 1992] — to the camera motion computation problem relies on 
interval arithmetic to take into account multiple constraints and screen space prop- 
erties. Satisfying camera movement parameters with respect to some given scene 
description — are obtained through constraint solving. The constraints involved are 
handled with Alg. J LA (Tab. II, p. 9). 

The main objective that motivated our work was to build a high-level tool allow- 
ing an artist to specify the desired camera movements for a "shot" using cinematic 
primitives. The resulting description is then translated into a constraint system in 
terms of the camera parameters, and solved using local consistency-based pruning 
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Fig. 4. Camera model 

techniques. A huge set of solutions is usually output as a result, from which a 
challenging task is to extract a limited representative sample for presentation to 
the user. 

The mathematical model used for representing the camera, its motion, and the 
objects composing a scene is presented in the following section; then follows the de- 
scription of some of the cinematic primitives together with their translation in terms 
of constraint systems. Addressing the problem of the isolation of a representative 
solution sample is deferred until Section 7.2. 

6.1 Modelling the Camera and the Objects 

A camera produces a 2D image by using a projection transformation of a 3D space 
scene. In the following, the image is referred to as the screen space (or image space) 
and the scene filmed as the scene space or 3D space. 

The standard camera model is based on Euler angles to specify its location and 
orientation. The work presented here is not bound to this representation, though, 
and any other representation would be convenient as well. 

A camera (Figure 4) possesses seven degrees of freedom, viz. its Cartesian position 
in the space, its orientation, and its focal length: 

— Position. Three scalars: x, y, and z; 
— View direction. Three scalars: 

—Pan. e, 

— Tilt, (p, 

— Roll. ip; 

— Focal length. One scalar: 7. 

Most movies are made of a large number of short elementary "shots" . Therefore, 
a simple model for camera motion may be adopted without losing too much expres- 
siveness. Sophisticated camera movements are obtained by assembling sequences 
of shots (Refer to Christie and al. 's work [Christie et al. 2002] for a way to perform 
this task). Due to a lack of space, the following description of camera movements 
(Figure 5) is restricted to primitive ones, and the reader is referred to the excellent 
book by Arijon [Arijon 1976] for a thorough presentation of the others: 
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Horizontal panoramic Travelling (tracking shot or doUy) Zoom-in and Zoom-out 



Fig. 5. Camera movements 

— Panoramic shot. A panoramic shot may be horizontal or vertical (i.e. around 
vertical or horizontal axis). The camera location is usually constant; 

— Travelling (tracking shot or dolly). A general term for a camera translation; 

— Zoom in or zoom out. A variation of the focal length of the camera. 

To our knowledge, existing declarative camera movement generators compute 
camera animation frame by frame [Druckcr 1994] or use calculated key frames [Shoe- 
make 1985] (fixed camera location and orientation), and interpolation of the camera 
parameters in-between. 

In contrast, the tool described here is based on a parametric representation. Com- 
puting a satisfying solution boils down to determining all the camera parameters 
such that some constraints on the scene are satisfied. Using the previous remark 
concerning simple elementary shots, parameters are modelled with degree three 
polynomials whose unknown is time t. For example, an horizontal pan (i.e. a move- 
ment along the horizontal orientation angle 6) might be defined as: 9{t) = cot + dg, 
where dg represents the initial horizontal orientation of the camera (at time t = 0) 
and eg is a constant velocity. Consequently, given Ve the maximum allowed pan 
velocity, one may note that cg must lie in the domain [—Vg, Vg], and then dg must 
lie in the domain [0, 27r], in order to make any orientation starting point eligible. 

In our view, the scene is considered as a problem data. Hence, from this point 
onward, every object composing a scene is assumed to have its location, orientation 
and movement already set by the user. 

Object properties rely on bounding volumes (location) and object axis (orien- 
tation): each object is bounded by a volume called a bounding box. Compounds' 
bounding box is the smallest box containing the botmding boxes of all the objects 
involved. In addition, bounding boxes may be associated to any set of objects 
(like a group of characters). Many geometric modellers provide such a hierarchy of 
bounding volumes. 

The orientation of any object on the image is determined by three vectors: Front, 
Up and Right (Figure 6). 

6.2 Properties and Constraints 

The translation of declarative descriptions of scenes into constraints is now pre- 
sented. Three kinds of desired properties may be distinguished: 

(1) Properties on the camera; 
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Fig. 6. Object axis and bounding volume 



(2) Properties on the object screen location; 

(3) Properties on the object screen orientation. 

Camera properties are a means to set constraints on the camera motion. Trans- 
lation of properties on the objects of a scene is described below. The notion of 
frame is first introduced as an aid to constrain an object location on the screen: a 
frame is a rectangle whose borders are parallel to the screen borders. A frame may 
be inside, outside, or partially inside the screen, though it is usiially fully contained 
in the screen. For special purposes, the frame size or/and location may be modified 
during the animation. 

With frames, an artist can define the precise projection zone of an object from 
the scene (3D) to the screen (2D). Our belief in a creation helping tool leads us to 
prefer this kind of soft constraint to the exact 2D screen location of the projection 
of a 3D point. 

The user defines frames (interactively or off-line), then chooses properties to 
apply to a frame and an object. For example, the statement (Figure 7): 

"The sphere s with center {xs{t),yfi{t), z^{t)) and radius r must be fully 
included in the frame f defined by the bottom-left point {xj{t), yj{t)) and 
the top-right point {x'^ (t) , y^ (t)) during the 20 seconds of the shot filmed 
by the camera c located at {xc{t),yc{t), Zc{t)) with orientation 6c{t) and 



is translated into the nonlinear constraint system of equations and inequations: 




with: 



(x{t) 

yit) 



'x,{t) - x,{t)) sm9,{t) + {y,{t) - y,{t)) cose,{t) 
Xs{t) — Xc{t)) cos9c{t) sin0c(t) + 

' {y,{t) - y,{t)) sin(/)c(t) sin0,{t) + {z^{t) - Zc{t)) cos(/)c(t) 
Xs{t) — Xc{t)) cos9c{t) cos4>c{t)+ 

' {y,{t) - j/c(i)) sineS) cos.^c(i) + {z,it) - Zc{t)) sin (t>c{t) 



< 



z{t) 
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Fig. 7. Projection in a frame 

where (^x{t),y{t), z{t)) is the projection of the center of s in the screen space at 
time t. 

Note that, in order to faU back to an instance of the guaranteed tuning problem, 
it is possible to discard the three equations by replacing x(t), y{t) and z{t) in the 
inequalities by the corresponding right-hand side. 

The constraints resulting from the translation of the declarative description of 
"shots" contain occurrences of the universally quantified time variable t. 

7. EXAMPLES AND BENCHMARKING 

A high-level declarative modeller tool for camera motion has been devised in order 
to validate the algorithms presented in Section 5. The prototype is written in C++ 
and Tcl/Tk; Figure 8 presents its graphical user interface: the animated scene to be 
filmed is displayed in a window together with the bounding volumes, while another 
window contains some previously drawn projection frames. The user constructs 
frames, selects objects, and assigns properties to objects in the scene. The output 
is a set of satisfying camera paths, and corresponding animations are shown in an 
output window. 

In the next section, we present some problems used to assess the quality of 
Alg. IPA^ (Tab. V) using Alg. IC02'^'^^' (Tab. IV and Def. 4.5) as the inner con- 
tracting operator parameter T (hereafter referred as IPABC). Comparisons with 
Alg. J LA (Tab. II) used by Jardillier and Languenou [Jardillier and Languenou 
1998] for the same kind of applications are produced and commented. Some tech- 
niques to speed-up the computation and improve the representativeness of the boxes 
output are also described. 

7.1 Description of the benchmarks 

Parabola; a curve fitting problem. This simple benchmark [Jardillier and 
Languenou 1998] corresponds to finding all the parabolas lying above a given line: 

yte[0,2]: at^ + bt + c^2t-l with o G [0, 1], 6 G [0, 1], c G [0, 1] 

Circle; a trivial collision problem. Benchmark Circle is a collision problem: 
given B a point moving along a circling path, find all points A such that the distance 
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Fig. 8. The declarative modeller tool 



between A and B is always greater than a given value. Benchmarks Circle2 and 
Circles are instances of the same problem with respectively 2 and 3 points moving 
round in circles. For only one circling point, we have: 

VtehTT,^]: [^£[-5,5] 

r- with \ y €z [—5, 5J 

^/ (ri sin t - a;)^ + (ri cost - yY ^ di \ di ^ 5 

where di is the mandatory minimal distance between A and B, and ri = 2.5 is 
the radius of i?'s circling path. 

Satellite; a collision problem. Given n satellites swivelling around a planet 
(Figure 9), we are looking for the parameters of an (n + l)th trajectory on which 
to put another satellite. Obviously, this trajectory must be such that the added 
satellite never collides with the already launched ones. 

The position of the ith satellite at time t is defined by: 

(Xi (t) \ / di cos 9t sin ujit + 
yi{t) = (sin T^i sin 6*,; sin (wjt + + cos cos (wit + (?!),;)) 
(^) / \ d-i (~ COS -04 sin 9i sin {tOit + 0; ) + sin V'i cos {uJit + (pi)'^ 

where variables 9i and (jji define the orientation of the plane of revolution, uji the 
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Axis of revolution 




Initial positu 
on the trajectory (b 



Fig. 9. A collision-free problem 



angular velocity of the satellite and (pi its initial position at t = 0. Variable di 
stands for the radius of the circle of revolution. Considering n satellites, wc arc 
looking for consistent values of the unknowns 9j , , iVj , of the satellite j = n+1 
such that: 



Vi G [-TT, tt] : < 



f dist(fi(t),fj(t) 
dist(f2(t),fj(t) 



> s 

> s 



dist(f„(t),fj(t)) ^ s 



where s represents the minimal distance allowed between two satellites. 

The unknowns to be computed are 9j , <f)j and ipj , with domains [0, 2'7r] . In this 
benchmark, we consider three satellites in the air with the following parameters: 



Parameter Satellite 1 Satellite 2 Satellite 3 



di 


5.0 


5.0 


5.0 


UJi 


1.0 


1.0 


1.0 


4>i 


0.0 


1.0 


2.0 


ei 


0.0 


1.0 


1.5 




0.0 


1.0 


1.5 



Robot; a collision problem. This benchmark is based on Example 1.1, p. 3. 
The actual constraint system to solve is: 



[0, 2] : ^ (a: - PAt)Y + {y - Py{t)y 



^ d 



where 

Pyit) 
ai{t) 
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di sinai(t) + d2 sin(ai(f) + a2{t) — tt) + da sin(ai(^) + a2{t) + 03(0) 
di cosai(t) + d2 cos(ai(t) + a2{t) - tt) + d3Cos(ai{t) + a2{t) + asW) 
t + w/A 
2t-l 
0.2t + 0.1 




The initial domains for the unknowns x and y arc both sot to [0,5], the size of 
the robot's hand is set to 0.5 and the respective lengths of the arm's segments are 
dx = 1.0, d2 = 2.0 and = 1.0. 

PointPath; a motion planning problem. This benchmark [Jaulin and Wal- 
ter 1996] introduces a simple motion planning problem (Figure 10). We need to 
compute an object's path, starting at position Mq and ending at position Mi, while 
avoiding collision with the ground and some objects. The path M{t) is tentatively 
chosen as a polynomial written as a linear combination of Bernstein polynomials of 
degree 3, and controlled by the points Pi and P2 of unknown coordinates: 

M(t) = MaBlit) + PiBl{t) + P2Bl{t) + MiBl{t) 

where the Bernstein polynomials arc: 

Blif) = (1 - f)\ Bf{t) = 3t{l i)\ Bl{i) = 'M\l - t), Bl = 

The constraints specify that for all t in the time frame: 

(1) M{t) must be above a curve representing the floor; 

(2) the distance between M{t) and a static object S = (4.8, 1.0) must be greater 
than 1. 

Which leads to: 

with M{t) = [ \. Domains for the coordinates of Pi and P2 are initialized to 



[-10,10]. 

Projection; a camera control problem. Benchmark Projection corresponds 
to the problem of projecting a moving sphere in a moving frame, already presented 
in Section 6.2. The initial values chosen for the camera are the following: Xc G 
[-3,3], 2/e G [-3,3], Zc = 2, G [-0.5,0.5], 9^ = 0, and 7e = 0.8. 

Projectioni is the original problem presented in Section 6.2; Projection^ is the 
same problem with one more constraint on the distance between the sphere and 
the camera; Projections, is the problem with two frames and two spheres. 
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GarlofFGrafl ; inner approximation computation [Garloff and Graf 
1999]. Find the values for the parameters v and w ensuring the stabiUty of the 
polynomial: 

p{s) = + vs'^ + {iv — 5v — 13)s + w 

Using the Linard-Chipart criterion, the problem can be transformed into the 
equivalent system: 

v,w > 

—5v'^ — 13v + vw — w > 

Following Garloff and Graf, we compute the inner approximation of the relation 
-5i)2 -13v + vw-w>0forv€ [2, 10] and w G [40, 50]. 

GarloffGraf2; inner approximation computation [Garloff and Graf 

1999]. This is again a stability problem: find sound values for A, B, and D veri- 
fying: 

A,B,D>0 

-AB + A + D'^-D-1>Q 
AB-AD-2A + D'^+ AD^ + 4D > 
AB^ - AB^D - 4AB^ + 2ABD + AAB + 2BD^ + 5BD^ + 2BD - D'' - AD^ - AD > Q 

AB-2A- BD^ - ABD - AB - 2D^ +3D-2>0 

Following Garloff and Graf, the initial domains chosen are A S [100, 120], B g [0, 2], 
and D G [10,20]. 

According to Abdallah and al. [Abdallah et al. 1996], a CAD-based software such 
as QEPCAD (http: //www. cs.usna.edu/~qepcad/B/QEPCAD. html) requires more 
than two hours to prove the existence of a solution to the problem. 

7.2 Improving Computation 

Solvers such as Numerica usually isolate solutions with variable domains around 
10~^ or 10~^^ in width. By contrast, the applications this paper focuses on are less 
demanding since the resulting variable domains are used in the context of a display 
screen, a "low resolution" device. In practice, one can consider that a reasonable 
threshold e for the splitting process during the search is some value lower or equal 
to 10-2 or 10-3. 

One of the drawbacks of Algorithm J LA [Jardillier and Languenou 1998] is that 
successive output solutions are very similar, while it is of importance to be able 
to provide the user with a representative sample of solutions as soon as possible. 
Tackling this problem using Alg. IPABC is done as follows: given a constraint 
system of the form Vt G It : ci A • • • A Cm and a Cartesian product of domains 
B = /i X • • • X J„, we have two degrees of freedom during tlic^ solving process, viz. 
the selection of the next constraint to consider, and the selection of the next variable 
to split. Figure 11 presents the differences with regard to the order of generation 
of solutions for Circle2 for two strategies concerning the variable splitting order: 
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Fig. 11. Depth-first vs. semi-depth-first 



— depth-first, where each constraint is considered in turn, and each domain is spHt 
to the threshold spUtting Umit; 

— semi- depth- first where each constraint is considered in turn, but each variable is 
split only once and then put at the end of the domain queue. 

As one may see, the semi-depth-first algorithm computes consecutive solutions 
spread over all the search space, while the depth-first algorithm computes solutions 
downward and from right to left. 

Some strategies on constraint consideration order have also been investigated, 
whose impact on speed is described hereunder. Four strategies may be singled out: 

Simple method. Box consistency is applied on the negation of each constraint in 
turn; 

Pre-parse method. This method interleaves the simple method with a pre-parse 
algorithm: given a constraint c, and t the universally quantified variable, k canonical 
intervals are extracted from the domain It of i, and the consistency of c is tested for 
every one of them. At this stage, a failed check is suSicient to initiate a backtrack; 

Normal method. Box consistency is computed both for each constraint and for 

its negation; 

Global method. Given constraints Ci,Ci-^-l, . . . ,Cm, the global method applies the 
normal method on Cj, then checks whether the output boxes are consistent with the 
remaining constraints by mere evaluation. 

Charts in Figure 12 present the time spent for obtaining the first and all solutions 
for four benchmarks on a SUN UltraSparc 1/167 MHz under Solaris 2.5. 

Considering first Chart A, one may see that the simple method is the most in- 
teresting strategy for computing the first solution, while the normal method is very 
time consuming for problems with "many" constraints (projections). On the other 
hand, while the pre-parse method is a bad choice for computing only one solution, 
it is competitive for obtaining all solutions. 
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Normai Simple Global Pre-parse A. Normal Simple Global Pre-parse B. 

Fig. 12. A) First solution. B) All solutions (times in seconds) 




Fig. 13. Impact of the time splitting threshold oj on precision in J LA (from left to right: ui = 0.5, 
w = 0.1, and u) = 0.05) 



7.3 Results 

Algorithms JLA and IPABC provide diflFerent sets of solutions for the same problem. 
Consequently, a direct comparison of their performances is quite difficult. Moreover, 
the actual implementation of JLA uses a splitting threshold lu for slicing the domain 
of the universally quantified variable t instead of checking consistency by eventually 
reaching canonicity of the samples of the domain It. Figure 13 shows the impact 
of the threshold on the computed solutions for benchmark Circle2. 

Table VI (resp. VII) compares algorithms JLA and IPABC from the speed point 
of view for computing the first solution (resp. all solutions). Times are given in 
seconds on a Pentium III at 800 MHz under Linux. 

As one can see, the efficiency of IPABC compared to the one of JLA increases 
steadily with the precision required. 

We were not able to obtain any result for PointPath with JLA after several hours. 
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Table VI. JLA vs. IPABC — First solution 



Bciiclimaik JLA IPABC JLA/IPABC 



L cUL diUUid 1 c — -LU J 




0.04 


0.02 


2.0 


Parabola (e = IQ-^) 




0.77 


0.02 


38.5 


Circle2 (e = IQ-^) 




3.52 


0.01 


352000 


Circle2 (e = 10"^) 




3.57 


0.01 


357000 


Satellite (e = lO^^) 




0.91 


0.99 


0.91 


Satellite (e = 10"^) 




68.7 


0.99 


68.48 


Robot (e = 10-1) 




0.13 


0.01 


13 


Robot (e = 10-2) 




0.50 


0.01 


50 


Projections (e = 10- 




0.05 


0.07 


0.71 


Projections (e = 10- 




0.14 


0.07 


2.0 


Projectiong (e = 10- 




3.01 


0.07 


43.0 


Projection4 (e = 10^ 




0.11 


0.03 


3.66 


Projcction4 (e = 10^ 




1.11 


0.03 


37 


GarlofTGrafl (e = 10 




15.79 


0.01 


1579 


GarlofrGraf2 (s = 10 




5.58 


0.04 


139 


PointPath (e = 0.5) 




??? 


7.85 


??? 



Table VII. JLA vs. IPABC— All solutions 



Benchmark JLA IPABC JLA/IPABC 



Parabola 




10.65 


1.22 


8.72 


Circle2 




3.16 


0.15 


21.06 


Circles 




3.41 


0.55 


6.2 


Satellite 




>600 


54.91 


> 175 


Robot 




22.65 


1.97 


11.5 


Projections 




>600.00 


3.43 


>175.4 


Projections 




182.09 


16.01 


11.35 


GarlofTGrafl (e 


= 10-2) 


422.54 


1.49 


283.6 


GarlofrGraf2 (e 


= 10-2) 


29.54 


1.04 


28.40 


PointPath (e = 


0.5) 


??? 


534.14 


??? 



8. CONCLUSION 

Unlike the methods used to deal with universally quantified variables described 
in [Hong and Buchberger 1991], the algorithms presented in this paper are purely 
numerical ones (except for the negation of constraints). Since they rely on "tradi- 
tional" techniques used by most of the interval constraint-based solvers, they may 
benefit from the active researches led to speed up these tools. What is more, they 
are applicable to the large range of constraints for which an outer c;ontracting oper- 
ator may be devised. By contrast, CAD-based methods deal with polynomial con- 
straints only, as is the case with the method based on Bernstein expansion [Garloff 
and Graf 1999]. 

Despite the dramatic improvement of the new method described herein over 
the one given by Jardillier and Languenou, handling of complex scenes with many 
objects and a camera allowed to move along all its degrees of freedom in a reasonable 
time is beyond reach for the moment. Nevertheless, a comforting idea is that most 
of the traditional camera movements involve but few of the degrees of freedom, 
thereby reducing the number of variables to consider. 
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Following the work of Markov [Markov 1995] and Shary [Shary 1995], a direction 
for future research is to compare the use of Kaucher arithmetic [Kaucher 1980] 
to compute inner approximations of relations with the use of outer contracting 
operators. Their work is also related to the one by Armengol et al. [Armengol et al. 
1998] and Gardenes et al. [Gardenes and Trepat 1980; Gardenes and Mielgo 1986] on 
modal interval arithmetic. Yet, these approaches require algebraizing trigonometric 
constraints, an operation known to slow down computation [Pau and Schicho 2000]. 

APPENDIX 

Set of floating-point numbers 
SmaUest floating-point number greater than a 
Greatest floating-point number smaUer than a 
Set of closed intervals whose bounds are floating-point numbers 
Smallest box containing p: Outer(p) = P\{B e I" | p C B} 
Inner approximation of p: 
Inner(p) = {r g R" | Outer({r}) C p} 
Union of the boxes in S: 

HBi,... ,Bn}5 = {r e R" I 3j e {!,..., n} s.t. r e Bj} 
Relation associated to the constraint c: 
Pc = {{ri,...,rn) I c(ri,...,r„)} 

Relation associated to Vii G / : c 

"Negation" of the constraint c: (/ ^ is defined as / ^ 0) 
Domain of the variable x in the box B 
Set difl^erence of Bi and Ba as a set of boxes: 
^Bi H B25 = {r G Bi I Outer({r}) n B2 = 0} 
Box B where the interval I is replaced by the interval J 
Power set of the set <S 
Split in k boxes the box B 

Split in k boxes the box B (never splits the interval corresponding to the 
domain of v 

union of vectors of sets componentwise: 

[jjm,...,K)} = i[jui...,\jK) 
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